ON EXACT SOLUTIONS AND NUMERICS FOR 
COLD, SHALLOW, AND THERMOCOUPLED ICE SHEETS 

ED BUELER AND JED BROWN 

Abstract. This three section report can be regarded as an extended appendix to (Bueler, Brown & 
Lingle 2006) . First we give the detailed construction of an exact solution to a standard continuum 
model of a cold, shallow, and thermocoupled ice sheet. The construction is by calculation of 
compensatory accumulation and heat source functions which make a chosen pair of functions for 
thickness and temperature into exact solutions of the coupled system. The solution we construct 
here is "Test G" in (Bueler et al. 2006) and the steady state solution "Test F" is a special case. 
In the second section we give a reference C implementation of these exact solutions. In the last 
section we give an error analysis of a finite difference scheme for the temperature equation in the 
thermocoupled model. The error analysis gives three results, first the correct form of the Courant- 
Friedrichs-Lewy (CFL) condition for stability of the advection scheme, second an equation for error 
growth which contributes to understanding the famous "spokes" of (Payne et al. 2000), and third 
a convergence theorem under stringent fixed geometry and smoothness assumptions. 



1. Derivation of the exact solution 

1.1. Review: Equations of the continuum model. The flat bed, frozen base case of the cold 
shallow ice approximation is, for the purposes of this paper, taken to be the equations in the 
"Continuum Model" part of (Bueler et al. 2006); that paper is hereafter referred to as "BBL." The 
notation used in, the physics of, and the boundary conditions for the continuum model are all laid 
out in BBL. The equations are repeated here for convenience: 

dH 

mass-balance: — — = M — V • Q, (1) 

dt w 

dT k d 2 T dT 
temperature: — - = — 5 — U • VT — w— — h S. (2) 

at pc p oz z oz 

effective shear stress: {& X zi &yz) = —pg(H — z)S7H, (3) 

constitutive function: F(T,a) = Aexp (~^pj °' n ~ 1 > ( 4 ) 

horizontal velocity: U = -2pgVH f F(T,a)(H - () d(, (5) 

Jo 

f H 

map-plane flux: Q = / Udz, (6) 

Jo 

vertical velocity: w = — V • U d(, (7) 
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strain heating: £ = — (a xz , a yz ) ■ (8) 



The primary scalar unknowns are the thickness H, which equals the surface elevation in this flat 
bed case, and the absolute temperature T. Note a = {o\ z + cr^) 1 ^ 2 . As in (Bueler et al. 2006), to 
construct exact solutions we must allow an additional source term in the temperature equation: 

at oz pc p oz z 

1.2. Review: Specification of geometry and temperature. We repeat the basic specification 
of the exact solution from BBL. Suppose 

H(t,r)=H s (r)+4>(r) 1 (t) (10) 

where H s (r) is given by 



H, 



H s(r) = , t - "° ?T ,^ (1 + l/n)s - l/n + (1 - s)^ n - s^ n " ' , (11) 



(1 - l/ n )«/(2n+2) 

and f(r), g(t) are given by 



l+l/n _ 1+1/n 



n/(2n+2) 



if 0.3L < r < 0.9L and </>(r) = otherwise. Also 

j(t) = A p s'm(2irt/t p ). (13) 

We suppose from here on that the size of the perturbation f(r)g(t) is limited so that the slope 
dH/dr in the positive radial direction is negative. 
Now suppose 

T(t , r ,^ T , (r) K^M 

where v is found from H and T s by 



, s kTJr) / / H(t,r)G\ , . 



Then T satisfies the boundary conditions 



T| 2 _ H = T s (r) and ^ = -G/A; 

2=0 



The form of T in equation (|14j) is primarily constrained by the need to do the integral for U 
analytically (below). 

1.3. Computation of the velocity field. The content presented so far has appeared in (Bueler 
et al. 2006), but now we start the more detailed computation. 
The horizontal velocity is found from equations (3), (4), and (5): 

U(t, r, z) = -2{pgTA\VH\^VH j" exp (^y) (H - () n d( 

= -f2( w ) n A| J ff r | n " 1 i/ r e-^ f e-**(H - () n d( 

Jo 

= -v1{pg) n A\H r \ n - l H r e- Q ^ BT ^ l r^ n+ ^ / e°d n dB, (16) 



ON EXACT SOLUTIONS AND NUMERICS FOR THERMOCOUPLED ICE SHEETS 



3 



where H r = dH/dr and after a change of variables in the integral. We have generally suppressed 
dependence on t and r in these and remaining equations; dependence on z or £ will remain explicit. 
Here we have also introduced 

^' r) := rt s (v + hy 

The change of variables mentioned above is 9 := fj,(h — Q. It allows us to rewrite 

rz rflH 

/ e~ K {H - C) n d( = /x- (n+1) e-^ / e e 6 n d6. (17) 

JO Ju(H-z) 



Note U =0. 

I 2 = 



Suppose n = 1, 2, 3, . . . is an integer. Then we can do integral in ()16|) analytically. In fact, let 
p n be the polynomial defined by the relation 



e x x n dx = p n {x) e x + c. 

It is easy to see by integration-by-parts that p n (x) = x n — np n ^\{x) and po 0*0 = 1. In particular, it 
follows that ps (x) = x 3 — 3x 2 + 6x — 6 andp^x) = x 4 — 4x 3 + 12x 2 — 24x + 24. Alsop n (0) = (— l) n n\. 
Assuming that H r < for all r, define 

io{t,r) := 2{pg) n A{-H r ) n e- Q / RT °fi- (n+ V 

and 

I n (t,r,z) := / e e B n d9. (18) 

Then we see, by the definition of p n , that 

I n (t, r, z) = Pn (nH)e» H - p n (fi(H - z))e^ H ^ 

and it follows that 

U(t,r,z) = +TuI n (z). 

We now have an analytical expression for U involving no integrals. 

Regarding the above calculation, we believe that for ice sheet flow the integer cases n = 1,2,3,4 
completely suffice for numerical testing as the range 1.8 < n < 4 is the broadest range of exponents 
in the constitutive relation known to the authors (Goldsby & Kohlstedt 2001). It may eventually 
be appropriate to consider noninteger n cases, in which case the "incomplete gamma function" 
enters (Abramowitz & Stegun 1965), but we do not see the need presently. 

Next we seek analytical expressions for the vertical velocity w and for the term V • Q which 
appears in equation (1). Both of these quantities are vertical integrals of the horizontal divergence 
of the horizontal velocity U. 

Recall that the divergence of a radial function is V • (f(r)r) = r~ 1 (d/dr) [rf(r)]. Also, from now 
on we will liberally use "/ r " as an abbreviation for df /dr. 
Thus 



Id. _ . 1 du dl n 

V • U = -— [rujl n \ = -UJl n + —I n + U)-^-. 

r or r or or 



But 



du 
dr 



H rr QT S , 1 n fJ>r 



noting that T s is a function of r only, and 
dT 
dr 
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For the last calculation, recall that if F(x) = f^r) <p(t) dt then F'(x) = g'(x)4>(g(x)) — f'(x)<p(f(x)). 
As mentioned, we will integrate V • U vertically. In order to facilitate this integration, define 

<Mt,r):=i+n|^ + |||-(n + l)^ and 7 (t, r) := fe^i^H + pH r )H n 

so that 

V • U = lo [<t>I n (z) - {ifef^nr) e~» z (H - z) n+l - {ii n+1 e» H H r ) e~^ z {H - z) n + 7] . 
It follows that 

w(t,r,z) = - I V ■\]dQ = uj[ii- 1 {ii- 1 ^-4>)ln^i{z) + {4>{H-z) + H r )I n {z)- 1 z\ . (19) 
Jo 

We have again used the change of variable <|17JI . Also, we have integrated I n by changing order of 
integration: 

pz f/J-H r/iH rz 

I n (()d(= / / e e 6 n d9d(= / e d 9 n ( / <%) d9 



Jn(H-t) Jfi(H-z) ^J(H-0/n) 

(Z - H)I n {z) + ^In+liz). 



Note w\ „ = 0. 

I z=0 



1.4. Computation of the compensatory accumulation and heating term. Next, 

V • Q = VH ■ U| + I V • UdC = Vfl" • U| _„ -w\ 



from equation (6) and by the above-mentioned rule for differentiating integrals. (This is the 
calculation which shows the equivalence of the vertically-integrated equation of continuity (1) and 
the surface kinematical equation (10).) Thus 

V Q = -oj^- 1 (/i~V - <f>) I n +i(H) + u-fH. (20) 

We can now compute the compensatory accumulation by using equations (|2()|) and QlOjl: 

u-°JL + V.Q. 

Next we get to the strain heating term. From equations (3), (4), and (5), 

Thus, from equations (3) and (8), 

E _ — { - M H - ,)V*) • £ - ^^exp (- -|f±^) - . 

pc p oz Cp V RT s (v + H)/ 

The above completes much of the hard work. Recalling equations (1) and (17), we find the 
desired compensatory heat source: 

dT „ m dT k d 2 T „ 
E c = -E7 + U • VT + w- — S. 
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It remains to use the chosen functions H, T in QlOj). (|T1|) . respectively, and find the derivatives of 
H and T which appear in the relevant PDEs. In particular, 

H t = fg', 

, (is t + H t )(v + z) - (y + R)v t 



{v + z 



{2 



_ ^ v + H , rr {vr + H r ){y + z)-(u + H)u r 



v + z ' (y + z)^ 

2T S 



(v + z) 3 ' 

where u Ff," u F r ," U F Z " U F ZZ " denote partial derivatives with respect to the given variables. (Recall 
that, however, "T s " denotes the surface boundary value for temperature. Note also that /, g, T s 
are functions of a single variable.) 

From the above analysis, we see that one must compute at least the following list of analytical 
derivatives: H' s , H'J, f, f", g' , T' s , v t , v T fj, r ,T t ,T r ,T z ,T zz . This is actually done in the 

reference implementation in the next section. 

2. Reference implementations of Tests F and G 

This section contains a C code which accepts t, r, z and computes H, M, T, U, w, S, and S c for 
the exact solution in each of Tests F and G. Both these C codes (and corresponding Fortran 90 
codes) are in the public domain under the GNU Public License and are available at 

www . dms . uaf . edu/ ~bueler/ icef lowpage . htm 

The code here is the authoritative, detailed description of the exact solution. In particular, it 
includes the constants used in BBL, and has analytically-expanded forms for all of the derivatives 
listed above. These codes have been compiled using the GNU gcc compiler. 

The code is not particularly written for efficiency and can undoubtedly be modified for speed, 
though perhaps at some loss in clarity. An efficiency is included, namely, the code allows a one- 
dimensional array of vertical coordinates z as input and returns corresponding arrays for all of the 
z-dependent quantities. One may therefore ask the subroutine for values of T or E c at every depth 
in a particular column of ice, for instance. A very simple example program, which merely evaluates 
the exact solution and prints the result to the standard output, is included below. 

2.1. exactTestsFG.c. Here is the actual code to compute the exact solutions. 



#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 



double p3 (double x) { 

/* p_3=x"3-3*x"2+6*x-6, using Horner's */ 
return -6.0 + x*(6.0 + x*(-3.0 + x)); 

} 



double p4 (double x) { 

/* p_4=x"4-4*x~3+12*x~2-24*x+24, using Horner's */ 
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return 24.0 + x*(-24.0 + x*(12.0 + x*(-4.0 + x))); 

} 

int bothexact (double t, double r, double *z, int Hz, 

double Cp, double *H, double *M, double *TT, double *U, 
double *w, double *Sig, double *Sigc) { 

/* 

int bothexact (const double t, const double r, const double z[], const int Mz, 
const double Cp, double &H, double &M, double TT[], double U[] , 
double w[], double Sig[], double Sigc[]) { 

*/ 

const double pi = 3.14159265358979; 

const double SperA=31556926 . ; /* seconds per year; 365.2422 days */ 

/* parameters describing extent of sheet : */ 

const double H0=3000.0; /* m */ 

const double L=750000.0; /* m */ 

/* period of perturbation; inactive in Test F: */ 

const double Tp=2000.0*SperA; /* s */ 

/* fundamental physical constants */ 

const double g=9.81; /* m/s~2; accel of gravity */ 

const double Rgas=8.314; /* J/(mol K) */ 

/* ice properties; parameters which appear in constitutive relation: */ 

const double rho=910.0; /* kg/m~3; density */ 

const double k=2.1; /* J/m K s; thermal conductivity */ 

const double cpheat=2009 . 0; /* J/kg K; specific heat capacity */ 

const double n=3.0; /* Glen exponent */ 

/* next two are EISMINT II values; Paterson-Budd for T<263 */ 

const double A=3.615E-13; /* Pa~-3 s"-l */ 

const double Q=6.0E4; /* J/mol */ 

/* EISMINT II temperature boundary condition (Experiment F) : */ 
const double Ggeo=0.042; /* J/m"2 s; geo. heat flux */ 
const double ST=1.67E-5; /* K m~-l */ 
const double Tmin=223.15; /* K */ 

const double Kcond=k/ (rho*cpheat) ; /* constant in temp eqn */ 

/* declare all temporary quantities; computed in blocks below */ 
double power, Hconst, s, lamhat, f, goft, Ts, nusqrt, nu; 
double lamhatr, fr, Hr, mu, surfArr, Uconst , omega; 
double Sigmu, lamhatrr, frr, Hrr, Tsr, nur, mur, phi, gam; 
double I4H, divQ, Ht , nut; 
double I4,dTt,Tr,Tz,Tzz; 
int i ; 
double *I3; 

13 = (double *) malloc(Mz * sizeof (double) ) ; /* need temporary array */ 
if (13 == NULL) { 

fprintf (stderr , "bothexact () : couldn't allocate memory\n"); 

return -9999; 

} 

if ( (r<=0) I I (r>=L) ) { 

printf ("\nERR0R: code and derivation assume 0<r<L !\n\n"); 
return -9999; 

} 
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/* compute H from analytical steady state Hs (Test D) plus perturbation */ 

power = n/(2*n+2) ; 

Hconst = H0/pow(l-l/n, power) ; 

s = r/L; 

lamhat = (l+l/n)*s - (1/n) + pow(l-s , 1+1/n) - pow(s , 1+1/n) ; 
if ((r>0.3*L) kk (r<0.9*L)) 

f = pow( cos(pi*(r-0.6*L)/(0.6*L)) ,2.0); 
else 

f = 0.0; 
goft = Cp*sin(2.0*pi*t/Tp) ; 
*H = Hconst*pow(lamhat , power) + gof t*f ; 

/* compute TT = temperature */ 
Ts = Tmin+ST*r; 

nusqrt = sqrt( 1 + (4 . 0* (*H) *Ggeo) / (k*Ts) ); 
nu = ( k*Ts/(2.0*Ggeo) )*( 1 + nusqrt ); 
for (i=0; KMz ; i++) 

TT [i] = Ts * (nu+(*H)) / (nu+z[i]); 

/* compute surface slope and horizontal velocity */ 
lamhatr = ( (1+1/n) /L) * ( 1 - pow(l-s,l/n) - pow(s,l/n) ); 
if ( (r>0.3*L) kk (r<0.9*L) ) 

fr = -(pi/(0.6*L)) * sin(2.0*pi*(r-0.6*L)/(0.6*L)) ; 
else 

fr = 0.0; 

Hr = Hconst * power * pow(lamhat,power-l) * lamhatr + goft*fr; /* chain rule */ 
if ( Hr>0 ) { 

printf ("\nERR0R: assumes H_r negative for all 0<r<L !\n"); 
return 1 ; 

} 

mu = Q/(Rgas*Ts*(nu+(*H))) ; 
surfArr = exp(-Q/ (Rgas*Ts) ) ; 
Uconst = 2.0 * pow(rho*g,n) * A; 

omega = Uconst * pow(-Hr,n) * surfArr * pow(mu,-n-l) ; 
for (i=0; KMz ; i++) { 

I3[i] = p3(mu*(*H)) * exp(mu*(*H)) - p3(mu* ( (*H)-z [i] ) ) * exp(mu* ( (*H) -z [i] ) ) ; 

U[i] = omega * 13 [i] ; 

} 

/* compute strain heating */ 
for (i=0; KMz ; i++) { 

Sigmu = -(Q*(nu+z[i] )) / (Rgas*Ts* (nu+(*H) ) ) ; 

Sig[i] = (Uconst*g/cpheat) * exp(Sigmu) * pow( fabs(Hr)*( (*H) -z[i]) ,n+l) ; 

} 

/* compute vertical velocity */ 

lamhatrr = ((1+1/n) / (n*L*L) ) * ( pow(l-s, (l/n)-l) - pow(s, (l/n)-l) ); 
if ( (r>0.3*L) kk (r<0.9*L) ) 

frr = -(2.0*pi*pi/(0.36*L*L)) * cos (2 .0*pi* (r-0 . 6*L) /(0 . 6*L) ) ; 
else 

frr =0.0; 

Hrr = Hconst*power* (power-1) *pow(lamhat ,power-2 . 0) * pow(lamhatr,2.0) + 

Hconst*power*pow(lamhat , power-1) *lamhatrr + goft*frr; 
Tsr = ST; 

nur = (k*Tsr/(2.0*Ggeo)) * (1 + nusqrt) + 

(1/Ts) * (Hr*Ts-(*H)*Tsr) / nusqrt; 
mur = ( -Q/(Rgas*Ts*Ts*pow(nu+(*H) ,2.0)) ) * ( Tsr*(nu+(*H))+Ts*(nur+Hr) ); 
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phi = 1/r + n*Hrr/Hr + Q*Tsr/(Rgas*Ts*Ts) - (n+1) *mur/mu; /* division by r */ 
gam = pow(mu,n) * exp(mu*(*H)) * (mur*(*H)+mu*Hr) * pow((*H),n); 
for (i=0; i<Mz; i++) { 

14 = p4(mu*(*H)) * exp(mu*(*H)) - p4(mu* ( (*H)-z [i] ) ) * exp(mu* ( (*H) -z [i] ) ) ; 

w[i] = omega * ((mur/mu - phi)*I4/mu + (phi* ( (*H)-z [i] )+Hr) *I3 [i] - gam*z [i] ) ; 

} 

/* compute compensatory accumulation M */ 
I4H = p4(mu*(*H)) * exp(mu*(*H)) - 24.0; 

divQ = - omega * (mur/mu - phi) * I4H / mu + omega * gam * (*H) ; 
Ht = (Cp*2.0*pi/Tp) * cos(2.0*pi*t/Tp) * f; 
*M = Ht + divQ; 

/* compute compensatory heating */ 

nut = Ht/nusqrt ; 

for (i=0; i<Mz; i++) { 

dTt = Ts * ((nut+Ht)*(nu+z[i])-(nu+(*H))*nut) * pow(nu+z [i] , -2 . 0) ; 

Tr = Tsr*(nu+(*H))/(nu+z[i]) 

+ Ts * ((nur+Hr)*(nu+z[i])-(nu+(*H))*nur) * pow(nu+z [i] , -2 . 0) ; 

Tz = -Ts * (nu+(*H)) * pow(nu+z [i] , -2 . 0) ; 

Tzz = 2.0 * Ts * (nu+(*H)) * pow(nu+z [i] , -3 . 0) ; 

Sigc[i] = dTt + U[i]*Tr + w[i]*Tz - Kcond*Tzz - Sig[i]; 

} 

free(I3) ; 
return 0; 

} 



2.2. simpleFG.c. This is the simple program to exercise the above code. 



#include <stdio.h> 
#include <stdlib.h> 
#include "exactTestsFG.h" 

int main() { 

const double SperA=31556926 . ; // seconds per year; 365.2422 days 

const double Cp=200.0; // m; magnitude of the perturbation in test G 

double year, r, HF, MF, HG, MG; 

double *z, *TF, *UF, *wF, *SigF, *SigcF, *TG, *UG, *wG, *SigG, *SigcG; 
double *mb; /* a block of memory */ 
int j , Mz ; 

printf ("Enter t and r separated by newline"); 
printf(" (in yrs and km, resp.; e.g. 500 500) :\n"); 
scanf (""/.If ",&year); 
scanf ("'/.If ",&r) ; 

printf ("Enter z values sep by newline (in m);"); 
printf (" >-l> to end; e.g. 100 500 1500 -l:\n"); 

z = (double *) malloc(501 * sizeof (double) ) ; 
if (z == NULL) { 

f printf (stderr, "simpleFG: couldn't allocate memory\n"); return -9999; } 



j=0; 
do { 
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scanf ("'/.lf",&z[j]); 

if (j>490) printf ( " \n\n\nWARNING : enter -1 to stop soon!!!\n"); 
} while (z[j-l]>=0.0) ; 
Mz=j-1; 

mb = (double *) malloc(10 * Mz * sizeof (double) ) ; 
if (mb == NULL) { 

f printf (stderr, "simpleFG: couldn't allocate memory \n" ) ; return -9999; } 
TF=mb; UF=mb+Mz*sizeof (double) ; wF=mb+2*Mz*sizeof (double) ; 
SigF=mb+3*Mz*sizeof (double) ; SigcF=mb+4*Mz*sizeof (double) ; 
TG=mb+5*Mz*sizeof (double) ; UG=mb+6*Mz*sizeof (double) ; 
wG=mb+7*Mz*sizeof (double) ; 

SigG=mb+8*Mz*sizeof (double) ; SigcG=mb+9*Mz*sizeof (double) ; 
/* evaluate tests F and G */ 

bothexact(0.0,r*1000.0,z,Mz,0.0,&HF,&MF,TF,UF,wF,SigF,SigcF) ; 
bothexact (year *SperA , r * 1000 . , z , Mz , Cp , &HG , &MG , TG , UG , wG , SigG , SigcG) ; 

printf ("\nResults:\n Test F Test G\n"); 

printf (" (functions of r (resp. t and r) only):\n"); 

printf (" H = %12.6f (m) H = '/,12.6f (m)\n" ,HF,HG) ; 

printf (" M = %12.6f (m/a) M = 7.12. 6f (m/a)\n", 

MF*SperA,MG*SperA) ; 
for (j=0; j<Mz; j++) { 

printf ( " (z='/.10 . 3f ) : \n" , z [ j] ) ; 

printf(" T ='/.12.6f (K) T ='/.12.6f (K) \n" , TF [ j ] , TG [ j ] ) ; 

printf(" U = , /.12.6f (m/a) U = '/.12.6f (m/a)\n" ,UF[j] *SperA, 

UG[j]*SperA) ; 

printf (" w = , /,12.6f (m/a) w = '/.12.6f (m/a) \n" , wF [j] *SperA, 

wG[j] *SperA) ; 

printf (" Sig = y.l2.6f (*) Sig = '/.12.6f (*)\n", 

SigF [j] *SperA*1000 . , SigG [j ] *SperA*1000 . 0) ; 
printf (" Sigc = , /.12.6f (*) Sigc = '/.12.6f (*)\n", 

SigcF [j ] *SperA*1000 . , SigcG [j ] *SperA*1000 . 0) ; 

} 

printf ("(units: (*) = 10"-3 K/a)\n"); 

free(mb) ; 
return 0; 

} 



A run of simpleFG looks like this: 

$ simpleFG 

Enter t and r separated by newline (in yrs and km, resp.; e.g. 500 500): 

500 

500 

Enter z values sep by newline (in m) ; to end; e.g. 100 500 1500 -1: 



100 
500 
1500 
-1 



Results : 
(functions 



Test F Test G 

of r (resp. t and r) only): 
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H 


1925 


295290 


(m) 


H 


2101 


899734 


(m) 


M 


-0 


010510 


(m/a) 


M 





040738 


(m/a) 


0.000) 
















T 


265 


122620 


(K) 


T 


267 


835036 


(K) 


U 





000000 


(m/a) 


U 





000000 


(m/a) 


w = 





000000 


(m/a) 


w = 





000000 


(m/a) 


Sig = 





264346 


(*) 


Sig = 


1 


215392 


(*) 


Sigc = 


-0 


373726 


(*) 


Sigc = 


-1 


323664 


(*) 


100.000) 
















T 


263 


137595 


(K) 


T 


265 


849860 


(K) 


u 





661716 


(m/a) 


U 


2 


244496 


(m/a) 


w = 





000005 


(m/a) 


w = 


-0 


000758 


(m/a) 


Sig = 





173915 


(*) 


Sig = 





817817 


(*) 


Sigc = 


-0 


306255 


(*) 


Sigc = 


-1 


022931 


(*) 


500.000) 
















T 


255 


486095 


(K) 


T 


258 


194962 


(K) 


U 


1 


785938 


(m/a) 


U 


6 


217140 


(m/a) 


w = 





000291 


(m/a) 


w = 


-0 


011984 


(m/a) 


Sig = 





028439 


(*) 


Sig = 





149934 


(*) 


Sigc = 


-0 


199905 


(*) 


Sigc = 


-0 


340039 


(*) 


1500.000) 
















T 


238 


172200 


(K) 


T 


240 


856843 


(K) 


u 


2 


036372 


(m/a) 


U 


7 


227603 


(m/a) 


w = 





002288 


(m/a) 


w = 


-0 


050018 


(m/a) 


Sig = 





000029 


(*) 


Sig = 





000400 


(*) 


Sigc = 


-0 


193301 


(*) 


Sigc = 





365908 


(*) 



(units: (*) = 10~-3 K/a) 

These numbers allow an easy check on correctness if modifications are made to the implementation 
exact solutions or, for example, upon recompilation on a new machine. The numbers can be 
generally compared to the figures in (Bueler et al. 2006). 

3. Stability and convergence of a numerical scheme for temperature 

3.1. A traditional finite difference error analysis. In this section we analyze the error in 
a semi-implicit, first-order-upwinded finite difference scheme for equation @, namely the scheme 
used to generate the results in BBL. Appendix A of that paper provides a description of the coupled 
numerical scheme which solves equations (^Q) and @. This section "fleshes out" Appendix B of 
that paper, which sketches the error analysis here. 

The most important caveats about the analysis here is that the components of the velocity field 
are assumed to be known functions independent of the temperature T and also that the geometry of 
the ice sheet is assumed fixed. In these ways we are not analyzing the coupled numerical scheme. 
Nonetheless we believe this analysis provides enough information to help build a reliable adap- 
tive time-stepping scheme and also reveals a significant point of error growth in thermocoupled 
circumstances. 

One can compare the material here to (Calvo, Diaz & Vazquez 2002) which does a finite element 
analysis of a moving geometry and velocity field problem for the temperature equation in a flow 
line model but for which the thermomechanics are only semi-coupled to the geometry. 

We generalize equations © and © slightly. In particular we analyze the equation 

T t + u(x, y, z, t)T x + v(x, y, z, t)T y + w(x, y, z, t)T z = KT ZZ + f(x, y, z, t, T), (21) 

for absolute temperature T. We denote the exact solution to (|21j) by "T" or ll T(x,y,z,t)." The 
function / in (|2*Tj) generalizes £ and £ + £c which appear in @ and @, respectively. Note that £ 
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depends upon the temperature, the effective shear stress, and the strain rates. Thus £ is a function 
of variables x, y, z, t, T because the velocity components are assumed to be known functions. 

We suppose (|21[) applies on some bounded region C I 3 which is fixed in time, and we as- 
sume a rectangular computational domain enclosing f2. Consider a regular, rectangular grid on 
that computational domain, in four variables (x,y,z,t), with spacing Ax, Ay, Az, At and grid 
points denoted (xi,yj, z k ,ti). Let UijM = u(xi,yj, Zk,ti), etc, and assume T-j k is our numerical 
approximation of T(xi ,yj,z k ,ti). 



The numerical scheme is 



T^-Tl jk UP (Tj jk u ijkl ) Up [Tl k v ijkl ) Up [Tl 



+ 



ijt 

+ 



Wijkl 



At Ax Ay Az 

T , ii,k+l ijk ' ij,k—l ,■ i , rj-,1 \ 

= K — ^2 — + f(xi,yj,z k ,ti,T ijk ) 



(22) 



where 



Up (ip. 



a. 



Oi((pi - (fi-i), at > 
ai(ip i+ i - ifi), on < 0. 

That is, the advection terms are upwinded in the easiest first-order manner and the vertical conduc- 
tion term is treated implicitly in the easiest centered-difference manner. Note that we abbreviate 
the term f(xi,yj,z k ,t h T^ k ) by in what follows. 

The first kind of error we analyze is local truncation error, but only because it plays a supporting 
role in the analysis of the approximation error (total numerical error). The local truncation error 
is defined to be the nonzero result of applying the finite difference scheme to the exact solution 
(Morton &; Mayers 2005). In our case the local truncation error depends upon a long list of cases 
for the upwind scheme, but listing these turns out to be unnecessary (and certainly uninteresting). 
We assume for the next few equations that Uijki, Vijki,Wijki are all nonnegative, and we will soon be 
able to return to expressions which apply to all upwinding cases. (Equation (1221) can be rewritten 
without the "Up(- •)" notation if the signs of the velocity coefficients are known, of course.) 

Denote the local truncation error at a grid point by Tj,-jy, so 

_ T(xj,yj, z k ,U +1 ) - T(xj,yj, z k ,U) T(xj,yj, z k ,U) - T(xj^i,yj,z k ,ti) 

TW— Af + u ijkl Ax 

T(xi,yj,z k ,ti) -T(xi,yj-i,z k ,ti) T(x if yj, z k ,ti) - T(xi,yj, z k - X ,ti) 

+ VijH 1 r Wijkl- 



K 



Ay 0J ™ Az 

T(xi,yj,z k+1 ,ti +1 ) - 2T(xi,yj,z k ,ti +1 ) +T(x i ,yj,z k ^ 1 ,ti +1 ) 



Az 2 

- f(xi,yj,z k ,ti,T(xi,yj,z k ,ti)). 

The finite difference quotients here all have well-known Taylor expansions. Because T(x, y, z, t) 
solves (|21[). it follows that 

T m = 0(At,Ax,Ay,Az). 
Including the Taylor expansions and all cases of upwinding, 

r m = ±T tt At ± ^fT,.,.A.r ± ^T yy Ay ± ^fr zz Az - ^T zzzz Az 2 (23) 

where the higher partial derivatives of T are evaluated somewhere in the grid neighborhood of 
(xi, yj, z k ,ti). We see that the finite difference scheme is consistent (Morton & Mayers 2005). Note 
that if Wij k i = then Ty k i = 0(Az 2 ), but this case is too special to be of interest. 
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Let 

e ijk = Tijk — T(xi, l/j, Zk,tl) 

be the (signed) error at a grid point. Because Tl- k solves difference scheme ((221) exactly, and by 
the definition of the local truncation error, 

^ijk ijk c ijk c i-l,jk c ijk c i,j-l,k c ijk c ij,k-l (0A s 
+ Uijkl T VVijkl T \-Wijkl T (24) 



At 1 >JM Ax ' lJM Ay ' lJM Az 

T~ijkl ■ 



J+l _ n J+l , J+l 
j ^ij,k+ l ^ijk ~ ij,k~l 



Az 2 



f(T(xi,yj,z k ,ti)) - f{T l ijk ) 



We are interested in the evolution of error so we solve for the error at grid point (xi,yj,Zk) at 
time ti + \\ 

KAt\ KAt ( \ f At At At] , , , 

1 + 2 A^ J = ( e Sfe+i + e !S-i) + I 1 " «««^ " Vi i kl ~Ay ~ Wi ^A~ z j (25) 



At j At f At ; 



At 



f(T(xi,yj,z k ,ti)) - f{T l ijk 



A-tTijM 



Note that the quantity in curly braces actually varies depending on the upwinding do all 

coefficients of the errors which depend on the velocities. The possibilities for the quantity in curly 
braces are described by 

At At At 

1 ± Uijkl^— ± %-fcZ-r— ± Wijkl-r- , (26) 
Ax Ay Az 

with, for example, a "— " in front of Uijki if Uijki > and a "+" if Uijki < 0. 

The next step is significant. We identify an assumption which is sufficient, under the noted 
additional assumptions of smooth and fixed velocity fields and geometry, but which is presumably 
not necessary, for the stability of our finite difference scheme. This assumption is part of a maximum 
principle argument for a finite difference scheme; see standard examples of such arguments in 
(Morton & Mayers 2005). 

Assumption 1. The space-time grid satisfies 

. At . At , . At , . 

1 - \ U ijkl\-^: - \ V ijkl\-^ ~ \ W ijkl\-^ > (27) 
or, equivalently, the time step is chosen small enough to satisfy 

At< (VhM \viM \vhM\~\ (28) 

~ V Ax Ay Az J v ' 

As is standard in these maximum principle arguments, the plan now is to use this assumption 
to conclude that all coefficients of errors are nonnegative in (|25jl. and thus to bound the error. Let 

E l = max |e' ?fc | and t\ = max |Tyfcj|. 
i,j,k J i,j,k 

In particular, E is the maximum absolute approximation error over the whole grid at time t/. It 
would be nice to know that it does not grow too large; such is our goal. Equation (j2Sj) and all of 
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its upwinding variations collectively imply 

if AA, 74-1, KAt f , , At , At . At i 

1 + 2 A^J * 2 A^ E + I 1 " ^A~x ~ j A~y ~ ^A-J E (29) 



At , At , At 



+ At 



Atr z . 



f{T{xu Vj,z k , ti)) - f{Tl jk ) 

The reader should note that we have used Assumption ^ to derive (|29|). 

We have not yet arrived at the desired error inequality. In fact, the size of the " \f(T(xi, yj, Zfc,t/))— 
f{T-j k )\ n term is very important, and we must make another assumption which amounts to knowing 
that the derivative of / with respect to the temperature T is not too big. Though this assumption 
is fully justified in the case of the strain-heating term £ for ice, the Lipshitz coefficient which 
appears (below) is interesting in connection to the "spokes" . 

Assumption 2. For the source term f in (|21[). there exists a bounded nonnegative function 
L(x, y, z,t) > such that 

\f(x,y,z,t,Ti) - f(x,y,z,t,T 2 )\ < L(x,y, z,t)\Ti - T 2 \. 

In particular, if the partial derivative df/dT exists and is bounded, and if we define 



L f (x,y,z,t)= max 



— {x,y,z,t,T) 



then we may take L = Lf by the Mean Value Theorem. Here l T m i n " and "T max " are lower and 
upper bounds, respectively, on the exact solution. 

Note that an obvious lower bound on the exact solution T is the minimum value of the surface 
temperature in standard ice sheet circumstances. An upper bound is the pressure-melting temper- 
ature, in practice. The actual polythermal nature of ice sheets is obviously ignored here; compare 
(Greve 1997). 

One may call L in Assumption [21 a local Lipshitz constant for / (as a function of T). Let 
Lijki = L(xi,yj, Zk,U) be its grid value. Let L = sup( xy z \ e ^ t>0 L(x,y,z,t) be the bound on L, 
the global Lipshitz constant. 

We apply Assumption [21 to Q29JI an d collect terms by noting that several on the right have 
coefficients which add to one. Finally we recall the definition of e\- k . The result is: 

(l + \e^\ < 2^£ m + E l + AtL ijkl \e\ jk \ + Atr ; . (30) 

This is a significant error inequality. 

3.2. A convergence theorem. One consequence of inequality (|30j) is the following convergence 
theorem. 

Theorem 1. Assume that the velocity components u,v,w are smooth. Assume that the geometry 
is fixed and is sufficiently smooth (so that the local truncation error Tij k i is indeed O(At)). If 
AssumptionsU\and\Q apply then the error grows at an (at most) exponential rate in the total model 
run time tf, times an 0(At) factor: 

E l <t f exp(t f L) f max r k ^j < t f exp(t f L) O(At). 
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Therefore the scheme converges as the time step At goes to zero. 
Proof. By taking maximums over the grid, inequality (|30|) implies 

E l+1 < (1 + AtL) E l + Atr h 
We initialize the scheme with the correct initial values so E° = 0. By induction, therefore, 

E l < At \ i (l + AtL) l T + (l + AtL) l ~ 1 r 1 + --- + (1 + AtL^n^ + n. 
Recall that (1 + x/n) n < expx for x > and n > 0. It follows that 



E l < At[ max Tf. 

k 



1 + AtL) 1 + fl + AtL) 1 1 + h f 1 + AtL) 1 + 1 



iv 



< At ^maxT fc ^ / (1 + AtL) 1 < lAt fmax|r|fcj ^1 + 

< I At ^maxrfcj exp(tjZ) < t/ exp(tjL) ^maxr^ 

as claimed. (Note that the total time tf is the number of time steps multiplied by the step size: 
tf = NAt. In particular, for each step I we have I At < NAt = tf.) □ 

This theorem is desirable but it is not directly useful. It is not surprising given the strong 
smoothness assumptions and it is not practical because the convergence constant "iyexp(ijL)" is 
almost certainly too large in practice. 

Dismissing this theorem for now, we see three intermediate results of importance, namely the 
two Assumptions, which arise reasonably naturally, and the inequality (|30|) . 



3.3. Adaptive time-stepping. Note that inequality (|28[) is used in the adaptive time-stepping 
numerical scheme described in BBL. By default the mass-balance and temperature steps are syn- 
chronized. In this synchronized case the stability condition for the mass-balance scheme — see 
BBL — and constraint (|28j) are combined to give the time step 

a. • j nio n f 1 1 V 1 ( r, V 1 (\ u ijkl\ . \Vijkl\ , KifcilVU 

At = mm < 0.12 x 2 -— r + maxDm , max — h . + . > 

1 \Ax 2 Ay 2 J V y V \ Ax Ay Az J J 

where Diji is the computed value of the diffusivity of the mass-balance equation (BBL). The special 
value "0.12" is essentially empirical, i.e. it results from testing in a variety of circumstances. The 
"max"s which appear in the above equation are computed over the spatial grid at time t[. 

For fixed velocity and diffusivity fields, and if Ax = Ay, we see that At = 0{Ax 2 ) as Ax -> 
because of the mass-balance condition. As the grid is refined it might be expected that the mass- 
balance constraint will inevitably become active, but in fact the maximum magnitude of Wijki 
near the margin causes constraint (j28|) to be active in many cases including some EISMINT II 
simulations and when using a digital elevation map for the surface elevation of Antarctica. (See 
BBL on the source of large vertical velocities.) For the better behaved exact Tests F and G, and 
for realistic ice sheet simulations wherein the geometry has been smoothed by evolution and the 
grid is significantly refined, we do indeed see the mass-balance condition most active. 

One can reasonably consider a scheme which tolerates violations by a modest factor of the CFL 
constraint (|28j) at a few grid points. The cost would be possible loss of accuracy at locations 
where the computation is not likely to be accurate anyway (e.g. at margins, near mountain ranges 
within realistic ice sheets, or at transitions from frozen to sliding base if that transition is too 
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abrupt anyway). The benefit would be that the time step for the whole sheet could be longer and 
computation times could be reduced. 

3.4. Emergence of spokes. Let us reconsider error inequality Q3UJI . Suppose that we look at the 
worst case location on the grid at time that is, at i, j, k such that \e^\ = E l+1 . In that case 
(HI says 

\ e ijk\ < E l + AtLijki \e\jk\ + Air; (applies at worst case grid point). (31) 

This inequality says that the error at the worst-case grid point can be no larger than the worst 
error over the grid at the previous time step plus two terms which are proportional to the time 
step. It is a statement of limited growth of error unless these two terms happen to be big at time 

The first of these terms depends on the local size of the strain heating term, or, more precisely, 
its variability with respect to temperature. The second, involving the maximum local truncation 
error over the grid, depends on the smoothness of the exact solution, or, more precisely, on the 
degree to which it does not satisfy the finite difference scheme. We believe that because of the free 
boundary nature of real ice sheet problems, the second truncation error term is in fact likely to be 
large at points. It causes only arithmetic growth in error, however, and it depends purely on the 
smoothness of the exact solution rather than nonlinear and nontrivially evolving errors on the grid. 
In any case we concentrate on the first term because we believe it is involved in the emergence of 
the spokes in the EISMINT II results. 

It is an easy calculation that if a flow law of form (J3J) is used then 

2AQ n+1 / Q \ 1 Q 

® exp I -— -2 = S "5^2 • ( 32 ) 



df 






dT 




dT 



pc p R \ RT J T 2 RT 2 



Thus we may use 

Q 

Lijkl — ^ijkl 



RT 2 
1X1 ijM 



for the local Lipshitz constant in (|31|) . 

If, in fact, the two-regime Paterson-Budd (1982) flow law is used then we need the larger activa- 
tion energy for warmer ice, Qr>263.i5 = 13.9 x 10 4 Jmol -1 , in our worst case estimate of the local 
Lipshitz constant, rather than the cold value Qr<263.i5 = 6.0 x 10 4 Jmol -1 . Thus the rate of error 
growth jumps up as one moves from relatively cold ice to warmer ice, precisely into the region of 
warm spokes in EISMINT II experiment F (Payne et al. 2000). (The Hooke (1981) relation should 
behave in roughly the same manner. Note that Payne and Bladwin (2000) produced spokes for 
EISMINT II experiment F using both the Paterson-Budd and Hooke relations, though the details 
of the spokes differ.) In any case the quantity \dH/dT\ has a characteristic spatial variation, illus- 
trated in BBL, which strongly suggests that locations where this quantity is large are locations of 
the emergence of (warm) spokes. 
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